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We study the nonlinear driven response and sliding friction behavior of the phase-field-crystal 
(PFC) model with pinning including both thermal fluctuations and inertial effects. The model 
provides a continuous description of adsorbed layers on a substrate under the action of an external 
driving force at finite temperatures, allowing for both elastic and plastic deformations. We derive 
general stochastic dynamical equations for the particle and momentum densities including both 
thermal fluctuations and inertial effects. The resulting coupled equations for the PFC model are 
studied numerically. At sufficiently low temperatures we find that the velocity response of an 
initially pinned commensurate layer shows hysteresis with dynamical melting and freezing transitions 
for increasing and decreasing applied forces at different critical values. The main features of the 
nonlinear response in the PFC model are similar to the results obtained previously with molecular 
dynamics simulations of particle models for adsorbed layers. 
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I. INTRODUCTION 

In recent years, considerable attention has been given 
to the study of driven adsorbed layers in relation to 
sliding friction phenomena between surfaces at the mi- 
croscopic levelir— . The nonlinear response of the ad- 
sorbed layer is a central problem for understanding ex- 
periments on sliding friction between two surfaces with 
a lubricant^ or between adsorbed layers and an oscillat- 
ing substrate^*-. Various elastic and particle models have 
been used to study the driven dynamical transitions and 
the sliding friction of adsorbed monolayers^— A fun- 
damental issue in modeling such systems is the origin 
of the hysteresis and the dynamical melting and freezing 
transitions associated with the different static and kinetic 
frictional forces when increasing the driving force from 
zero and decreasing from a large value, respectively. In 
driven lattice systems, hysteresis can occur in the under- 
damped regime where inertial effects are present or in the 
presence of topological defects such as dislocations and 
thermal fluctuations^. Topological defects can be auto- 
matically included in a full microscopic model involving 
interacting atoms in the presence of a substrate potential 
using realistic interaction potentials. However, numeri- 
cal observation of the full complexities of the phenomena 
is severely limited by the small system sizes that can be 
studied, even when simple Lennard- Jones potentials are 
employed 3 . 

Recently, a phase-field-crystal (PFC) model was 
introduced^— that allows for both elastic and plastic 
deformations within a continuous description of the par- 



ticle density while still retaining information on atomic 
length scales. By extending the PFC model to take into 
account the effect of an external pinning potential- 3 , 
a two-dimensional version of the model has been used 
to describe commensurate- incommensurate transitions in 
the presence of thermal fluctuations^ and the driven re- 
sponse of pinned lattice systems without thermal fluctu- 
ations and inertial effects^. However, in order to study 
fully the sliding friction behavior, both inertial effects and 
thermal fluctuations need to be taken into account. 



In this work, we study sliding friction of an adsorbed 
layer via the nonlinear response of the PFC model to an 
external force in the presence of a pinning potential. To 
this end, we first derive the general stochastic dynamical 
equations for the particle and momentum density fields 
taking full account of the thermal fluctuations and iner- 
tial effects. The resulting coupled equations are studied 
numerically. At low temperatures, we find that the veloc- 
ity response of an initially commensurate layer shows hys- 
teresis with dynamical melting and freezing transitions 
for increasing and decreasing applied forces at different 
critical values. The main features of the nonlinear re- 
sponse are similar to the results obtained previously with 
molecular dynamics simulations of particle models. How- 
ever, the dynamical melting and freezing mechanisms are 
significantly different. In the PFC model, nucleation oc- 
curs via stripes rather than closed domains found in the 
particles model. 
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II. DYNAMICS OF THE PFC MODEL WITH 
INERTIAL EFFECTS 

In the generalized PFC model, the system is repre- 
sented by a coarse grained effective Hamiltonian that is 
a functional of the number density field. To take into 
account inertial effects in the dynamics, we need to con- 
sider the contribution of the kinetic energy to the to- 
tal energy of the system in addition to the configura- 
tional energy. Thus we consider the momentum den- 
sity, g(x) = p(x)v{x), a dynamical variable in the coarse 
grained Hamiltonian in addition to the particle density 
field p(x). The total effective Hamiltonian in the pres- 
ence of an external force / can be written as 



H t = H kln + Hi n t - / dxp(x)x ■ /, 



(1) 

where H k in is the kinetic energy contribution given by 
H kln = I dxf^S, (2) 



J 2p(x) ' 

and Hi n t (p) is the configurational contribution to the ef- 
fective Hamiltonian in the original PFC model. The last 
term is due to the presence of the the external force /. 

In the absence of energy dissipation, the time depen- 
dence of p and g are determined by the Poisson brack- 
ets {H t ,p} and {H tl g}. At finite temperatures, addi- 
tional dissipative noise terms are present in the dynam- 
ical equations. The noise satisfies the fluctuation dis- 
sipation relation s 16 ' 17 which allows the system to reach 
thermal equilibrium in the absence of external perturba- 
tions. We include the dissipative noise term directly in 
the dynamical equations for g, which is a nonconserved 
field. 

For the particle density p(x, t) we have 



and for the momentum density g(x, t) 

~Y1 J dS '{9i{x),gj{x')} 



(3) 



SH t 

-rj— 1- Ui{x,t), 



Sgj(x') 



(4) 



Substituting Eqs. © into Eqs. © and (@J gives 
dp 



where 7/ is a dissipative coefficient and the noise v(x, t) 
has variance 

(u i (S,t)u j (x',t')) = 2k B Tr]5(x - x" )S(t - t')5 id . (5) 

The Poisson brackets for the mass and momentum den- 
sities are given b y 16 ' 17 

{p{x), 9l {x')} = ViipfflStf - &)); (6) 
{9i{x),p{x)} = p(x)\7iS(x - x); 

{ gi {x), gj (x')} = V 3 -(ft(2)5(£ - £')) -V % ( 9l {x)5{x - x')). 



dt 



(7) 



— = -pVi— — - V—+pfi+Vi(x,t) 8 
at op p 

_ v^vt 9i9j 

^ J p ' 

j 

for a spatially uniform external force (/,). We can rede- 
fine the coefficient r\ pr\ to remove the denominator 
from the term gi/p in Eq. (jHJ. With this change the 
variance of the noise v{x, t) becomes 

(v i (x > t)v j (x f ,t')) = 2k B Tr]p{x)S(x ~ x')S(t - t')6 id (9) 

To leading order in the momentum density <?, we drop 
the quadratic term in Eq. © giving 



dp _ 



dg_ 

dt 



-B = -pv. 



SH int 

~jp~ 



+ Pfi - V9i + Vi(x,t). 



(10) 



(11) 



Similar dynamical equations for PFC models with inter- 
nal dissipation were obtained in Ref. HTJ. If the effective 
Hamiltonian Hint is known, then these coupled stochas- 
tic dynamical equations should provide a full description 
of the particle and momentum densities in presence of 
fluctuations represented by the noise fi(x, t) with corre- 
lations proportional to the temperature T and inertial 
effects determined by the damping parameter rj. In the 
overdamped limit, dg/dt = 0, with T = and / = 0, the 
equation for the time evolution of p obtained by inserting 
g from Eq. (fTTj) in Eq. (fit))) reduces to the deterministic 
equation for the density which has been obtained from 
classical density functional theory of liquids 21 . 



III. PFC MODEL WITH PINNING AND 
THERMAL FLUCTUATIONS 

In the presence of an external pinning potential^—, 
a specific form of the configurational energy contribu- 
tion Hint to the total effective Hamiltonian in Eq.(p} 
has been proposed which is an extension of the stan- 
dard PFC model free energy functional used in many 
application a 11 ' 12 . In dimensionless form, this effective 
interaction Hamiltonian Hint = ^pfc can be written as 



H p{c = J rfaf{^[r+(l + V 2 



+ ^ + (12) 



where ip(x) is a continuous field, V(x) represents the ex- 
ternal pinning potential and r is a parameter. The phase 
field ip{x) can be regarded as a measure of deviations of 
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the particle number density p(x) from a uniform reference 
value poj such that ip(x) = (p(x) — Po)/po- It is a con- 
served field and its average value, ip, represents another 
parameter in the model. The intrinsic wave vector of the 
model ki has no preferred directions and its magnitude 
is set to unity in the present work. 

In the absence of a pinning potential, the Hamiltonian 
of Eq. (fT2|) is minimized by a configuration of the field 
ip(x) forming a hexagonal pattern of peaks with recipro- 
cal lattice vectors of magnitude \kh\ ~ 1, in an appro- 
priate range of values for the parameters r and ip in the 
model. This periodic structures of peaks in ip can be 
regarded as a simple model of an atomic layer. 

To study the nonlinear dynamical behavior we take 
the form of Hint given by Eq. (fT2l) together with the ki- 
netic energy and the external force terms for the total 
effective Hamiltonian H t in the dynamical equations Eq. 
(fTTj) and in Eq. ([9]). We make the additional simplify- 
ing approximation that p(x) w p in Eq. ^ and in the 
coefficient of the first term in Eq. (|lip . This approxima- 
tion ensures that the effective diffusion coefficient in the 
model is positive definite for any temperature and driv- 
ing force. Another motivation for this approximation is 
to show that the dynamical equations used in the previ- 
ous work s 13 ' 15 follows from the more general Eqs. (p~0|) 
and (Hip . Setting po = 1, we obtain, 

f~V. 9 ; (.3, 
% „ SH pfc 

~dt = ~5xb~ ~ v i\ x i *)> 

(v i (x,t)v j (x',t')) = 2k B Ti 1 5{x - x')S{t - t')Sij. 

Here we have redefined g — > g + pof /rj to remove a uni- 
form term on right hand side of the equation for g. 

The above coupled equations can be combined in a 
single equation by applying the operator V- to both sides 
of the second equation and using the first one to eliminate 
V • g, giving 

(£(2,t)£(&,lf)) = 2k B Tr 1 V 2 8{x~x')8{t-t'). 

When the driving force / and the external pinning po- 
tential are set to zero, the dynamical equation above is 
identical to the one used in Ref. 22 to study propagat- 
ing density modes in the PFC model. It can also be 
obtained by introducing inertial effects in the dynamical 
equation of the PFC model through a memory function 
of exponential form^ 3 -. In the limit of large rj when dgi/dt 
in Eqs. (Ti"3l) or, equivalently, d 2 ip/dt 2 in Eqs. ([13]). 
can be neglected, these equations reduce to the familiar 
overdamped dynamical equations used in the previous 
work o 13 ' 15 without inertial effects at zero temperature. 



IV. NUMERICAL RESULTS AND DISCUSSION 

In this section, we present our numerical results for 
the velocity response of the PFC model in presence of 
an external pinning potential under a uniform applied 
force. For the numerical calculations, the phase field ip(x) 
and momentum density field g(x) are defined on a space 
square grid with dx = dy = ir/4 with periodic boundary 
conditions. System sizes L x L with L = 64 to 128 where 
used. The Laplacians and gradients were evaluated using 
finite differences. 

We consider a pinning potential V{x) representing a 
substrate with square symmetry 

V(x) = -Vo[cos(fc x) + cos(fc J/)], (15) 

where fco defines the period of the pinning potential for 
both the x and y directions. The lattice misfit between 
the phase-field crystal and the pinning potential can be 
defined as S m = (1 — fco). We choose the parameters 
of the PFC model as r = -0.25, $ = -0.25, a lattice 
mismatch 6 = —0.5 and a pinning strength Vq — 0.15. 
For these parameters, the ground-state configuration in 
the absence of an external driving force is a c(2 x 2) 
phas e - 14 ' 15 , where every second site of the lattice of the 
pinning potential corresponds to a peak in the phase field 
ip(r). This commensurate configuration is stable below 
a commensurate melting temperature T c » 0.055. This 
is determined from the temperature dependence of the 
structure factor as well as that of the specific heat. The 
structure factor S(k) is calculated from the the positions 
Rj of the peaks in the phase field as 

N P 

S (£) = <E — e-*^-^')). (16) 
w'=i P 

Here, (...) denotes a time average which is equivalent to 
a thermal average at equilibrium. For increasing tem- 
peratures thermal fluctuations disorder the layer and the 
scaled structure factor S(k c )/N p evaluated at the pri- 
mary reciprocal lattice vector fc c of the commensurate 
phase decreases from a value of unity at T = rapidly 
through T c to zero at higher temperatures. The transi- 
tion is broadened due to finite size effects as shown in 
Fig. [IJa). Another signature of this commensurate melt- 
ing transition is that the specific heat develops a broad 
peak near T c corresponding to an increase in the fluctu- 
ations in -0 at the transition as shown in Fig. [ljb). Fi- 
nally, the mobility p defined as \mif^Q(v/ f) where v is 
the drift velocity, also changes qualitatively through the 
transition. In the commensurate phase, a finite threshold 
for the sliding of the overlayer exists and hence the mo- 
bility is vanishingly small. As the system melts above T c 
the mobility rapidly rises and reaches a plateau at higher 
temperatures. This behavior is shown in Fig. [T};. For the 
study of the nonlinear response and sliding friction of the 
system, we focus on an initial state which corresponds to 
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a well pinned c(2 x 2) phase initially, corresponding to 
T = 0.01 <C T c , and a damping parameter rj = 0.4. 

In our numerical calculations, a driving force along the 
x direction is increased from zero to a maximum value 
above a critical depinning force and then decreased back 
to zero. For each value of the force, the coupled equa- 
tions (|13p are solved using an Euler algorithm with time 
step dt = 0.001 — 0.005. The implicit trapezoidal method 
was also used to check numerical stability. For each value 
of the force, 10 6 time steps were used to allow the sys- 
tem to reach a velocity steady state and an additional 
period with the same number of time steps were used to 
calculate the average velocity and other time averaged 
quantities. 

To study the velocity response of the PFC model, we 
need to determine the velocity of the peaks^ in the phase 
field ip(x). This is done by determining the time depen- 
dence of the peak positions (Ri(t)) in ip. The steady state 
drift velocity v for the system is obtained from the peak 
velocities Vi (dRi/dt) as 



N P 



(17) 



where Np is the number of peaks and (...) denotes time 
average. The steady state structure factor S(k) which is 
a measure of translational order, is also calculated from 
the the peak positions Rj, as in Eq. fT^)) . 

The most notable qualitative feature of the velocity 
response to the driving force / (as shown in Fig. [2]) is 
that it shows hysteresis behavior with two different crit- 
ical force thresholds /„ « 0.075 for increasing forces and 
fb ~ 0.045 for decreasing forces. These two threshold val- 
ues correspond to the static frictional force and kinetic 
frictional force respectively. As the force is increased be- 
yond f a , the velocity jumps abruptly from zero to a finite 
value whereas when the force is decreased below /& the 
velocity of the sliding layer drops abruptly and become 
pinned by the external potential to form an immobile 
commensurate state again. Below, we present the micro- 
scopic details of the configurations for the system in the 
neighborhood of the two thresholds. This allows us to 
characterize the change in velocity response at f a as a 
dynamical force induced melting transition of the initial 
commensurate state, and the second transition at f as 
a dynamical force induced freezing of the sliding phase. 

To study these transitions we first examine the behav- 
ior of the steady state structure factor, as shown in Fig. 
[3l We focus on the dependence of S(Q) on the driving 
force, where Q is the dominant reciprocal lattice vector 
of the layer. Q — k c corresponds to the primary recip- 
rocal lattice vector for the c(2 x 2) phase and Q = kh 
to the reciprocal lattice vector of the hexagonal phase in 
absence of the driving force. Consistent with the veloc- 
ity response behavior in Fig. [5J on increasing the force 
beyond f a , S(k c ) drops abruptly to zero. This is the 
onset of the dynamical force induced melting transition 
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FIG. 1: Temperature dependence of the (a) scaled structure- 
factor peak S(k c )/N p ; (b) specific heat C, and mobility jj, for 
the model without an external driving force. 
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FIG. 2: Velocity response as a function of applied force. Ar- 
rows correspond to critical values f a , fb and f c . 



of the initial c(2x2) commensurate state. The behavior 
of S(k c ) is analogous to the temperature induced disor- 
dering transition shown in Fig. [2(a). On decreasing the 
force, the value of S(k c ) stays vanishingly small until the 
force drops below the threshold fb, at which point S(k c ) 
rapidly increases to a value corresponding to the com- 
mensurate pinned state. The other interesting feature is 
that for / > f c w 0.12, the structure factor shows clear 
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FIG. 3: Scaled structure-factor peak S(Q) /N p as a function of 
applied force. Here Q stands for the primary reciprocal lattice 
vector for either the c(2 x 2) phase (k c ) or the hexagonal phase 
(kh)- Filled and open symbols correspond to increasing and 
decreasing forces, respectively. 
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FIG. 5: Fraction of density peaks pk with coordination num- 
ber k for decreasing applied force. 




6-fold coordinated peaks at the reciprocal lattice vectors 
Q = kh corresponding to a hexagonal phase, which grows 
as the driving force increases. This implies that at a 
driving force larger than this third threshold / c , there 
is another dynamic continuous transition from a disor- 
dered phase into an incommensurate hexagonal phase. 
This state emerges as the average effect of the external 
pinning potential becomes less and less important at high 
sliding velocities and the steady state then approximately 
corresponds to the phase-field crystal in the absence of 
the pinning potential which has hexagonal symmetry in 
the equilibrium state. However, since the scaled struc- 
ture factor for the hexagonal phase is still much less than 
unity, this incommensurate hexagonal phase is not fully 
ordered even at at the largest force values (/ < 0.15) 
studied so far. 
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FIG. 7: Configuration of the density peak locations with 
corresponding coordination numbers for the density plot in 
Fig. 
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To better understand the dynamic melting transition 
at f a and the freezing transition at we inspect the 
steady state configurations near the two thresholds as 
well as the region between the two thresholds in real 
space, which yield more direct and detailed information 
than the structure factor at the peak values. Treating 
the peaks in the phase field as " particles" , the coordina- 
tion number of each particle in the commensurate c(2 x 2) 
phase is 4 while for the incommensurate hexagonal phase 
it is 6. Furthermore, the core of dislocations in a hexago- 
nal lattice correspond to particles with 5-fold and 7-fold 
coordination numbers. Thus, further insight into the mi- 
croscopic nature of the steady state configurations can be 
obtained by inspecting the fraction of particles with co- 
ordination numbers pi, ps, ps, pi- The results are shown 
in Figs. U and Fig. [3] for increasing and decreasing ap- 
plied forces respectively. Fig. [4] shows that when the 
force / is increased beyond f a , besides the rapid drop in 
Pi consistent with the structure factor data, the fraction 
of other coordination numbers ^5, pe, pj also increases 
significantly untill / reaches f c , beyond which p 5 and p 7 
start to decrease and we have a continuous transition to 
an incommensurate hexagonal phase. Thus the nature 
of steady state above f a is a strongly disordered state 
analogous to the high temperature phase in the absence 
of driving force above the commensurate melting tem- 
perature. On decreasing the force below / Q , the data in 
Fig. \5\ shows that the system remains in a melted state 
with large disorder until the threshold fb is reached, be- 
low which we have only 4-fold coordination number and 
the system returns to a pinned c(2 x 2) phase. Taken 
together, the qualitative behavior of the structure factor 
and the coordination number strongly suggest that the 
transitions at f a and fb can be regarded as a force in- 
duced dynamical melting and freezing transition respec- 
tively. 

Finally, we look at a snapshot of the phase field i^(x) 
obtained in the steady state for a driving force just above 
the dynamical freezing threshold This is shown in 
Fig. O It consists of stripes of commensurate c(2 x 2) 
phase separated by disordered domain walls. These do- 
main walls are mobile liquid like regions as confirmed by 
a count of the fraction of various coordination numbers 
shown in Fig. [7J The commensurate stripes have only 4- 
fold coordination numbers, while in the liquid like domain 
walls, there is a mixture of 5-fold, 6-fold and 7-fold coor- 
dination numbers. Starting from this steady state which 
still has a non-zero average sliding velocity, we can watch 
the dynamical freezing in real time as the driving force is 
reduced to a value just below the freezing threshold fb- 
The time sequence of the freezing is shown in Fig. |8] In 
returning to the pinned state, the domain wall regions 
gradually shrink and eventually disappear. 

The main features of the dynamical melting and freez- 
ing and hysteresis effects in PFC model described above 
are similar to those found earlier by molecular dynamics 
simulations of particle models with interacting Lennard- 
Jones potentials 3 - and the uniaxial Frenkel-Kontorova 
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FIG. 8: Snapshots of the density field for increasing times 
after starting from the moving state of Fig. [6] and reaching 
the pinned commensurate state, at a force f x = 0.045 just 
below fb- Times are in units of 10 3 time steps. 

model^ under a driving force starting from a commensu- 
rate state. The similarity of results from these very differ- 
ent models demonstrates the universality of the hysteresis 
loop for models with inertia effects and the macroscopic 
consequence of stick and slip motion. When compared 
with the Lennard- Jones particle model one notable dif- 
ference is the mechanism of the dynamical freezing at fb- 
In the present PFC model, it involves parallel liquid-like 
domain walls similar to the uniaxial Frenkel-Kontorova 
model, although in the absence of the driving force there 
is no easy direction. In the Lennard- Jones model, nu- 
cleation and growth occur via closed pinned domains 3 . 
The origin of this intriguing difference requires further 
investigation of both atomistic and PFC models. One 
possibility is that the nucleation of stripes is related to 
the absence of a fixed constraint on the number of peaks 
in the PFC model. In this case, the mechanism of the dy- 
namical freezing found in the present PFC model should 
be compared to the results of particle models with a con- 
stant chemical potential rather than with a fixed particle 
number. Unfortunately, such results are currently un- 
available. 



V. CONCLUSIONS 

In this paper, we have derived general stochastic dy- 
namical equations for the particle and momentum den- 
sity fields including both thermal fluctuations and iner- 
tial effects. The new equations are applied to the study 
of the nonlinear response to an external driving force 
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for a PFC model with a pinning potential. The model 
describes a driven adsorbed layer as a continuous den- 
sity field, allowing for elastic and plastic deformations. 
The numerical results showed that at low temperatures, 
the velocity response of an initially commensurate layer 
shows hysteresis with dynamical melting and freezing 
transitions for increasing and decreasing applied forces 
at different critical values. The inclusion of both ther- 
mal fluctuations and inertial effects are crucial for a cor- 
rect description of these dynamical transitions. The main 
features of the nonlinear response, in particular the hys- 
teresis loop separating the static friction and sliding fric- 
tion thresholds are similar to the results obtained previ- 
ously with particle models. However, the details of the 
dynamical melting and freezing mechanisms are signifi- 
cantly different. In the PFC model considered here, they 
correspond to nucleation of stripes rather than closed do- 
mains found in particle models. It is possible to describe 
more realistic sliding adsorbed systems if the parameters 
of the model are adjusted to match experimental systems, 



similar to the recent works for the colloidal systems^ and 
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